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SUDDEN  STRETCHING  OF  A  FOUR- LAYERED  COMPOSITE  PLATE 


by 

6.  C.  Sih 

Institute  of  Fracture  and  Solid  Mechanics 
Lehigh  University 
Bethlehem,  Pennsylvania  18015 

and 

E.  P.  Chen 
Sandia  Laboratories 
Albuquerque,  New  Mexico  87115 

ABSTRACT 

A  research  effort  primarily  concerned  with  the  understanding  of  laminated 
composite  plates  with  cracks  subjected  to  time-dependent  extensional  loads  is 
reported  here.  When  loads  are  applied  suddenly  to  a  laminate,  waves  are  re¬ 
flected  and  refracted  through  the  laminae  and  give  rise  to  stresses  and  strains 
throughout  the  composite  system.  The  process  is  three-dimensional  in  character 
and  presents  a  formidable  problem  in  the  theory  of  elastodynamics,  particularly 
in  the  presence  of  crack-like  imperfections. 

An  approximate  theory  of  laminated  plates  is  developed  by  assuming  that  the 
extensional  and  thickness  mode  of  vibration  are  coupled.  The  mixed  boundary 
value  crack  problem  of  a  four-layered  composite  plate  is  solved.  Dynamic  stress 
intensity  factors  for  a  crack  subjected  to  suddenly  applied  stress  are  found  to 
vary  as  a  function  of  time  and  depend  on  the  material  properties  of  the  laminate. 
Stress  intensification  in  the  region  near  the  crack  front  can  be  reduced  by  having 
the  shear  modulus  of  the  inner  layers  to  be  larger  than  that  of  the  outer  layers. 
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INTRODUCTION 


The  current  interest  in  laminates  for  structural  application  is  associated 
with  the  high  strength-to-weight  ratio  which  can  be  developed  in  laminates. 

These  laminates  are  generally  composed  of  layers  which  have  been  reinforced  by 
embedding  unidirectional  fibers.  The  layers  are  adhered  to  each  other  such  that 
the  fiber  direction  varies  from  one  layer  to  the  next  in  a  previously  determined 
manner.  The  freedom  of  choice  for  fiber  orientation  in  the  layers  of  the  com¬ 
posite  system  enables  the  development  of  laminates  with  special  preferential 
directional  properties  for  particular  applications.  Because  of  this  character¬ 
istic  of  fibrous  composites,  the  employment  of  these  systems  rather  than  equiva¬ 
lent  homogeneous  members  will  be  clearly  advantageous  in  many  applications. 

Because  of  the  complicated  internal  structure  of  composite  systems,  stress 
analysis  is  much  more  difficult  than  for  equivalent  single-phase  material.  One 
fact  which  emerges  very  clearly  from  laminate  studies  is  that  the  stress  field 
in  composite  systems  is  truly  three-dimensional  in  character.  Thus,  even  the 
stress  field  in  a  symmetric  laminate  subjected  to  in-plane  loading  cannot  be  ac- 
curately  modeled  by  standard  two-dimensional  methods  of  analysis.  The  previous 
work  in  this  area  further  indicates  that  relatively  little  effort  has  been  made 
to  formulate  laminate  plate  theories  that  can  effectively  solve  for  the  redistribu¬ 
tion  of  stresses  and  strains  due  to  the  presence  of  mechanical  imperfections  such 
as  cracks. 

One  possible  means  of  simplifying  the  three-dimensional  equations  of  elas¬ 
ticity  is  to  invoke  the  concept  adopted  in  the  formulation  of  plate  theory.  Ap¬ 
proximate  stress  and  strain  dependence  on  the  plate  thickness  coordinate  are  as¬ 
sumed  such  that  the  governing  differential  equations  possess  only  two  independent 
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space  variables.  In  addition,  special  attention  must  be  given  to  the  state  of 
affairs  near  the  crack  when  formulating  plate  theories  for  analyzing  crack  prob¬ 
lems.  With  this  in  mind,  Hartranft  and  Sih  [1]  developed  an  approximate  three- 
dimensional  theory  for  a  single  material  plate  containing  a  through  crack.  The 
condition  of  plane  strain  was  preserved  ahead  of  the  crack  as  suggested  by  Sih 

[2] .  This  theory  was  later  extended  to  laminates  by  Badaliance,  Sih  and  Chen 

[3]  to  solve  the  problem  of  a  through  crack  in  a  laminar  plate  subjected  to  in¬ 
plane  loading.  The  through  crack  configuration  represents  a  preliminary  effort 
to  model  the  damage  of  composite  plates.  Additional  complications  arise  when 
the  load  is  time  dependent.  These  considerations  will  be  taken  into  account  in 
the  development  of  a  new  dynamic  theory  of  laminated  composite  plates  subjected 
to  extensional  loads. 

This  work  is  concerned  with  the  formulation  of  a  dynamic  theory  of  laminated 
plates  and  reduces  to  that  of  Kane  and  Mindlin  [4]  for  the  single  material  plate. 
The  idealized  condition  of  stress  and  displacement  continuity  across  the  inter¬ 
face  is  replaced  by  assigning  certain  conditions  of  material  nonhomogeneity  in 
the  thickness  direction  of  the  laminated  plate  as  if  it  were  a  single  layered 
nonhomogeneous  plate.  The  nonhomogeneity  is  made  equivalent  to  a  symmetric  lami¬ 
nate  balanced  with  reference  to  its  mid-plane.  A  through  crack  is  assumed  to 
exist  in  a  four-layered  laminate.  Dynamic  stress  intensity  factors  are  obtained 
for  the  case  of  a  suddenly  applied  uniform  in-plane  loading  and  shown  to  vary 
as  a  function  of  time.  Discussed  are  also  the  influence  of  material  properties 
of  the  layers  on  the  local  stresses. 
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BASIC  FORMULATION 


The  elastodynamic  equations  of  generalized  plane  stress  are  adequate  only 
if  the  frequency  of  vibration  is  lower  than  that  of  the  first  thickness  mode  and 
the  wave  length  is  large  in  comparison  with  the  plate  thickness.  In  other  words, 
the  coupling  between  extensional  and  thickness  mode  of  vibration  can  be  ne¬ 
glected.  When  laminated  composite  plates  are  stressed  dynamically,  loads  are 
transmitted  through  the  laminae  by  the  reflection  of  thickness  refraction  of 
stress  waves.  The  mode  of  vibration  cannot  be  ignored,  particularly  in  the  vi¬ 
cinity  of  a  crack-like  imperfection  where  the  stress  state  acquires  a  three-di¬ 
mensional  character. 

A  dynamic  laminate  plate  theory  will  be  developed  to  solve  the  problem  of  a 
four-layered  composite  plate  with  a  through  crack  subjected  to  a  suddenly  applied 
uniform  extensional  load.  The  theory  is  a  generalization  to  that  given  by  Kane 
and  Mindlin  [4]  for  a  single  layer  plate  in  which  the  extensional  and  thickness 
mode  of  vibration  are  assumed  to  be  coupled.  Accounted  for  is  the  lowest  thick¬ 
ness-stretch  mode  such  that  the  displacement  is  normal  to  the  plate  surface. 
Mindlin  and  Medick  [5]  have  also  considered  a  formulation  in  which  the  thickness- 
shear  mode  of  vibration  with  displacement  parallel  to  the  plate  surface  is  also 
included.  The  mid-plane  of  the  plate  is  taken  as  the  nodal  plane  of  vibration. 

Consider  a  four-layered  composite  plate  of  thickness  h  as  shown  in  Figure  1 
where  each  layer  has  the  same  thickness  h/4.  The  two  outer  layers  have  material 
properties  (y2>v2)  or  (X2>U2^  while  the  two  inner  layers  have  material  properties 
(u.|,vi)  or  (Xpu-]).  The  Lame  coefficients  are  denoted  by  Xj  and  (j  =  1,2). 

The  layers  are  stacked  such  that  symmetry  prevails  across  the  mid-plane  of  the 
laminate  composite.  The  crack  of  width  2a  cuts  through  the  entire  thickness  of 
the  laminate.  » 


X 
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Figure  1  -  A  four-layered  composite  plate  with  a  crack 


=  (X+2y)e^  +  X(ey+K:e2) 

<^y  =  (A+2p)ey  +  X(e^+K:e2) 

=  (X+2p)ic2e^  +  XK(e^+Ey) 
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Tyz  = 


■^ZX  = 


^xy  '^'^xy 


The  constant 


K  =  ir/  /iY 


accounts  for  the  coupling  between  the  extensional  and  thickness  mode  of  vibration. 

It  is  determined  from  the  three-dimensional  equations  of  elasticity.  As  in  the 

development  of  plate  theories,  the  resultant  strain  quantities  (a  )  ,  (A  )  , 

j  ^  j 

....  (A^^)  (j  =  1,2)  will  be  defined: 
xy 


[(A^)  ,  (Ay)  ,  (A^)  ,  (Ayy)  ]  =jf  /  [£^,Ey,€y,Yyy]dZ 


(''zy  ‘VV  '  h  'J4 


I 

+  /  [ey.Sw. 


X’  y’  z’ 'xy- 


"  "F  u-(z.  •^^xz’^yz^^'^^ 


Jin  ^^xz’^yz^^^^  .  •(,  ^^xz’^yz^^^^ 


Substituting  equations  (2)  into  (5),  it  is  found  that 


(Ax)^  -  (Aj^)^ 
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fj  V 

(Ay)^  -  (Ay)^  - 


9v  av 

(a  )  =  (A  )  =  —  +  — ^ 

'  '  xy''2  3y  9x 


2  ^''2 


2  ^''2 

(Ayz)^  -  (\z)^  ~  h  ^ 


The  laminate  plate  theory  can  be  most  conveniently  formulated  in  terms  of  the 
stress  resultants 


(N^.Ny.N^.N^y)  = 


and  the  transverse  shears 


(Rv»U  =  L  (^xz’''yz^^^^ 


The  stress-strain  relations  in  equations  (3)  when  enforced  yield 


L  o  V  U  V 

N^(x,y,t)  =  I  [(e+2Y)  ^  +  3  ^]  +  3kv. 


I  O  V  O  V 

Ny(x,y,t)  =  2  [(e+2Y)  +  3  j^h  +  Bkv^ 

,  9V  9V 

N2(x,y,t)  =  (3+2y)k^V2  +  J  3Kh 
,  3V  3V 

Nxy(x,y,t)  =  2  Yh  (-^  +  ■^) 
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1 

Rj^(x,y,t)  =  ^  — 


R„(x,y,t)  =  ttp  6h2 


48  3y 


In  equations  (9)  and  (10),  B,  y  and  6  stand  for 


3  =  X.j+X2»  Y  =  vii+P2»  =  yi+7vi2  (11) 

Denoting  and  p2  as  the  mass  density  of  the  inner  and  outer  layers  of  the  lami¬ 
nate,  the  three  equations  of  motion  governing  N  ,  N  , R  become 

^  y  y 

3x  3y  2  W~ 

3N  3N  ,  32v 

3R  3R  ,  32v_ 

^  3/  ■  "  48 


The  result  of  substituting  equations  (9)  and  (10)  into  (12)  renders 


„  3V„  3V,,  -JQ  3v,  32v^ 

(PlV  3t^ 


3V„  3V,  3^V„ 

YV^Vy  +  (3+y)  “F^  "  ^Pl'^P2^  ^ 


8V..  3V. 


^  (6.2y).^v,  -  =  (p,.7p3)  ^ 
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where  +  9V9y^  is  the  Laplacian  operator  in  two  dimensions. 


METHOD  OF  SOLUTION 


Equations  (13)  will  be  solved  by  introducing  two  potential  functions  (i)(x,y,t) 
and  H(x,y,t)  as 


Vx(x,y,t)  =  If  +  ly 


Vy(x.y.t)  =  If  -  |y 


Making  the  appropriate  algebraic  manipulations,  the  governing  equations  for  the 
potential  functions  can  be  derived  by  enforcing  equations  (13): 


=  (p^+Pj)  -  S(6+2t)7'*4i  +  ,2^ 


"  ((p,+P2)[(pi+7p2)  I#  +  ^  (B«y)  0] 


-  C5(Pi+p2)  +  (3+2y)(pi+7p2)  |p-  (v2((.)]} 


Once  (j)(x,y,t)  and  H(x,y,t)  are  known,  v  and  v  can  be  obtained  from  equations 

X  y 

(14)  and 


V2(x,y,t)  -  [(P1+P2)  lyl"  "  (B+2K)v2(j)] 


Suppose  that  a  uniform  stress  resultant  is  applied  suddenly  to  the  crack 
surfaces  and  the  resulting  deformation  is  symmetrical  with  respect  to  the  x-axis, 
then  the  following  conditions  are  to  be  specified: 
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(17) 


Ny(x,o,t)  =  -NpHCt),  x<a 

Vy(x,o,t)  =  0,  x^a 

where  H(t)  is  the  Heaviside  unit  step  function.  The  condition  of  symmetry  fur¬ 
ther  requires  that 

N^y(x,o,t)  =  Ry(x,o,t)  =  0,  for  all  x  (18) 

Use  will  now  be  made  of  the  Laplace  transform.  Let  ((>  (x,y,p),  H  (x,y,p),  etc., 
denote  the  Laplace  transforms  of  the  functions  (j)(x,y,t),  H(x,y,t),  etc.  Equa¬ 
tion  (15)  when  expressed  in  the  Laplace  transform  domain  become 

(v2-aj|)(|,*(x,y,p)  =  0 

(v2-a)2)(^*(x,y,p)  =  0  (19) 

(v2-(o2)H*(x,y,p)  =  0 

where  the  potential  (j){x,y,t)  has  been  separated  into  two  parts: 

<|)(x,y,t)  =  (j)^(x,y,t)  +  (|>2(x,y,t)  (20) 

in  terms  of  time  t  or 

ic  'ff  ' 

<!>  (x,y,p)  =  <J>i(x,y,p)  +  (|)2(x,y,p)  (21) 

in  terms  of  the  Laplace  transform  variable  p.  The  parameters  u.  (j  =  1,2,3)  in 

J 

equations  (19)  are  defined  as 
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1)?  o  =  [(oi  +6  )(-^)  +  p  ±  {[(a  +6  )(-^)  +  p  ] 

i,d  no  0  0  (0-  0  0  O'  OJ-  0 


n  2  n  2  V2 

-  4a  5  (M  [(■^)  +  p  ]}  1 

0  0 


‘>2  ~  (pi'*’P2)p^y 


in  which  the  newly  defined  quantities  are 


3+2  Y  f.  .  6 

!p1+P2)Yq’  0  "  (p^+7p2)Yq 


^  4(pi+P2)  2  ^  12  k^(b+2y1 


P1+7p2  ’  h^(p^+p2j 


and  Yq  takes  the  form 


2  =  . . iiLg+x).— 

'o  (P1+P2) (3+2Yj 


Equations  (19)  then  give 


yc  ✓ 

<l'l(x,y,p)  =  -  /  A(s,p)  cos(sx)  exp{-s^y)ds 


yc  V 

^2(x.y»p)  =  -  /  B(s,p)  cos{sx)  exp{-S2y)ds 


H  (x,y,p)  =  ^  j  C(s,p)  sin(sx)  exp(-Soy)ds 
0 


with  s.  being  given  by 

J 


Sj  =  /s^+uj4,  j  =  1,2,3 
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The  dynamic  problem  has  now  been  reduced  to  finding  the  three  unknown  functions 
A(s,p),  B(s,p)  and  C(s,p). 

DUAL  COUPLED  INTEGRAL  EQUATIONS 

Before  the  boundary  and  symmetry  conditions  can  be  enforced,  it  is  necessary 

'k  k 

to  obtain  v  {s,y,p),  v  (x,y,p),  etc.,  in  terms  of  the  unknowns  in  equations 
X  y 

(25).  With  the  help  of  equations  (14)  and  (16),  it  can  be  shown  that 


v*(x,y,p)  =  -  7  /  [sA(s,p)  exp(-s,y)  +  sB(s,p)  exp(-S2y) 

X  TT  0  ' 

+  S2C(s,p)  exp(-S2y)]  sin(sx)ds 


Vy(x,y,p)  =  -  f  /  [s^A(s,p)  exp(-s^y)  +  S2B(s,p)  exp(-S2y) 


(27) 


+  sC(s,p)  exp(-S2y)]  cos(sx)ds 


V2(x,y,p)  =  7  /  [a^A(s,p)  exp(-s^y)  +  A2B(s,p)  exp(-S2y)]  cos(sx)ds 
0 


The  quantities  A.  (j  =  1,2)  are  given  by 

J 


4  ,  h(B^2Y)  .  „  ,2 

2eK  3+2y  “j-*’  ^ 


(28) 


k  k 

Similarly,  the  Laplace  transform  of  N^(x,y,p),  N  (x,y,p)  become 

A  y 


*  “  (Pl+Po)P^ 

Nj^(x,y,p)  =  7  yh  /  <-1—^ - sf]  A(s,p)  exp(-s^y) 


(p-|'*'Po)P^ 

+  [  2  -  s|]  B(s,p)  exp(-S2y).-  SS3C(S,P)  exp(-S3y)}cos(sx)ds 
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•][A(s,p)  exp(-s^y)  +  B(s,p)  exp(-S2y)] 


*  0  “  (pi+Po)P^ 

Ny(x,y,p)  =  f  yh  /  {[s^  +  — ^ - 


+  ss^CCsjP)  exp(-S2y)}cos(sx)ds 


N*  (x,y,p)  *  I  yh  /  {ss^A(s,p)  exp(-s^y)  +  ss2B(s,p)  exp{-S2y) 


+  ^  (s2+s2)  C(s,p)  exp(-S3y)}sin(sx)ds 


while  R^(x,y,p)  and  R  {x,y,p)  take  the  forms 
X  y 


R^Cx.y.p)  =  -  IJ:;;:  /  s[A^A(s,p)  exp(-s^y)  +  A2B(s,p)  exp(-S2y)]  sin(sx)ds 

(30] 

R*(x,y,p)  =  -  I57/  [s^a^A(s,p)  exp(-s^y)  +  S2A2B(s,p)  exp(-S2y)]  cos(sx)ds 


The  symmetry  conditions  in  equations  (18)  when  applied  show  that  A(s,p),  B(s,p) 
and  C(s,p)  can  be  expressed  in  terms  of  a  single  unknown  D(s,p): 


A(s,p)  =  D(s,p) 

n 

(p-l+Pp)P^  (Pl+P2)p^  0-,  ,  X 

B(s,p)  =  -s-j  [  ^  3+2y  ”  (■- 

2ss,  (a)?-a)2) 

- 

Application  of  the  mixed  boundary  conditions  in  equations  (17)  leads  to  a  sys¬ 
tem  of  dual  integral  equations 
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/  D(s,p)  cos(sx)ds  =  0,  x^a 
0 


II  n 

/  sF{s,p)  D(s,p)  cos(sx)ds  =  -  p-r— ,  x<a 


The  function  F(s,p)  is  known: 


S2+S| 


(p1+p2)P^. 


Hs.p)  =-^[Cs2  {1  - 


^  ( P-j  '*'p2  i 


2  2 


^,(p1+P2)P" 

e+2Y 

77v^ 

^2^  3+2y 


3+2y)-(A)^ 


The  standard  procedure  by  Copson  [6]  may  be  applied  to  solve  equations  (32)  and 


the  result  is 


D(s,p)  =  - 


irNoa^  ^  [(Pl+P2)P^/(2'^2Y)]-a)^  ^ 

^  '(ei+PzIP"  V  ,,  , 

- y -  <’  ■  S+27''"r“2> 


X  /  /?  $  (5,p)  J  (sa5)dc 
0  ° 


in  which  $  (?>p)  can  be  computed  from  a  Fredholm  integral  equation  of  the  second 


$  (?,P)  +  /  <5  (n,p)  L(c,n,p)dn 
0 


The  kernel  L(^,n»p)  is 


L(c»n,p)  =  ^  /  s[G(|,  p)  -  1]  Jo(s5)  J  (sTi)ds 
0 
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while  the  function  6(s,p)  is  related  to  F(s,p)  in  equation  (33); 


{ [  ( P-i +P2 )  P^/ (  3+2y  )  ]-a)^  }y 


DYNAMIC  STRESS  INTENSITY  FACTORS 


Of  interest  is  the  intensification  of  the  dynamic  stresses  ahead  of  the  crack. 
Hence,  the  integrals  in  equations  (29)  and  (30)  must  be  evaluated  for  large  values 
of  s  which  corresponds  to  distances  near  the  crack  edge  x  =  ±a  and  y=0.  In  terms 
of  the  polar  coordinates  r  and  e  in  Figure  1,  the  Laplace  transform  of  the  stress 
resultants  for  small  r  are  found: 

* 

*  (p)  A  A  n 

N^(r,e,p)  = - cos  |-  (1  -  sin  sin  +  0(r°) 


^(p) 


,(r,e,p)  =  — - cos  |-  (1  +  sin  |-  sin  |^)  +  0(r°) 


"k 

It  (P)  Q  Q  Oq  rt 

N^,,(r,0,p)  =  — —  sin  cos  -s-  cos  #  +  0(r°) 
xy  ^  ^  ^ 


R*(r,e,p)  =  R*(r,e,p)  =  0(r°) 


yc 

in  which  ki(p)  is  the  Laplace  transform  of  k^(t): 


k*(p)  = 


The  Laplace  inversion  theorem  may  now  be  applied  to  give 
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Equations  (40)  reveal  that  dynamic  loading  does  not  affect  the  functional  re¬ 
lationship  of  r  and  9.  The  stress  intensity  factor,  however,  is  a  function  of 
time: 


k-jCt)  = 


Nov^ 

Zui 


I 

Br 


exp(pt)dp 


(41) 


where  Br  denotes  the  Bromwich  path  of  integration.  Once  $  (?,p)  is  calculated 
from  equation  (35)  and  evaluated  at  C=1 >  equation  (41)  may  be  solved  numerically. 

Figure  2  gives  a  plot  of  $  (l,p)  as  a  function  of  C2i/pa  where  C21  =  (ui/p-])  ' 
is  the  shear  wave  velocity  referred  to  the  material  in  the  inner  layers.  For 
Pi  =  P2,  ■^‘i  =  ’^2  ”  *^1  ~  ^2’  increase  monotonically 

with  C2-j/pa-  Three  different  ratios  of  ^2/1^1  “  considered. 

Making  use  of  the  results  in  Figure  2,  ki(t)  in  equation  (41)  may  be  computed. 
Refer  to  Figure  3  for  a  display  of  ki(t)/NQya  versus  C2it/a,  The  resultant 
stress  intensity  factors  are  observed  to  vary  as  a  function  of  time.  Their  am¬ 
plitude  rise  quickly  reaching  a  peak  and  then  declines.  The  solution  for  a  homo¬ 
geneous  plate  corresponds  to  y2/y]  ~  1*0  Poisson's  ratio  and  mass  density 
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=  5.0 


C21  t/a 

Figure  3  -  Normalized  resultant  stress  intensity  factor  versus  Cp,t/a  for 
a/h  =1.0 
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for  the  inner  and  outer  layers  are  assumed  to  be  equal.  The  peak  value  of  k^{t) 
is  greater  than  that  of  the  homogeneous  plate  solution  for  y2>iii  while  the  opposite 
is  found  for  Hence,  the  intensity  of  the  crack  border  stress  field  can 

be  reduced  by  having  the  shear  modulus  of  the  outer  layers  to  be  smaller  than 
that  of  the  inner  layers. 


CONCLUDING  REMARKS 

A  dynamic  laminate  plate  theory  has  been  developed  for  solving  crack  bound¬ 
ary  value  problems.  The  complexity  of  the  problem  owing  to  material  nonhomogeneity 
and  dynamic  stress  analysis  necessitates  certain  simplifying  assumptions  so  that 
effective  analytical  solutions  can  be  obtained.  It  is  shown  that  the  dynamic 
stresses  near  a  mechanical  imperfection  such  as  a  crack  are  intensified  depending 
on  the  stacking  sequence  of  the  laminae.  In  general,  this  intensity  tends  to  in¬ 
crease  quickly  for  small  time  reaching  a  peak  and  then  decreases  to  the  static 
solution  for  sufficiently  long  time.  When  the  modulus  of  the  outer  layers  are 
smaller  than  that  of  the  inner  layers,  the  crack  border  stress  intensity  reaches 
a  maximum  quicker  than  the  homogeneous  solution  but  with  a  smaller  magnitude. 

The  opposite  holds  for  the  case  when  the  outer  layers  are  stiffer  than  the  in¬ 
ner  layers.  Information  of  this  type  is  useful  for  evaluating  the  resistance  of 
laminate  plates  to  impact  loadings. 
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COMPUTER  PROGRAM:  DYNAMIC  LAMINATE  PLATE  THEORY  WITH  A  CRACK 


1  PROGRAM  FLAP (INPUT. OUTPUT) 

real  non (4) »F(4.4,?)fG(4*4)*0(4).PT(4) 
real  8(4)*C(4) 
real  LP(19) .OTA (19) 

5  EQUIVALENCE  (NON.R) 

COMMON  K1.K?.K3.K4 
COMMON/AUX/H.P.PKl ,PK2.BMU.X.Y 
LP(1)=0.0 
OTA(1)=0.0 

10  read  2.K1.K2.K3.K4 

2  FORMAT (12)  . 

*  K1  =  ORDER  OF  SYSTEM  OF  EQUATIONS 

*  K2  -  NO.  OF  DISTINCT  KERNELS 

*  K3  =  NO.  OF  DATA  POINTS 

15  *  K4  *  NO.  OF  DATA  SETS  TO  BE  EVALUATED 

*  SET  UP  DATA  POINTS 

AK=K3 

DO  5  N=1.K3 
ANsN 

20  5  PT(N)=AN/AK 

*  SET  UP  INTEGRATION  MATRIX 

MSK3-2 

N=K3-1 

A=K3 

25  A=1,/(3.*A) 

DO  10  K=2.M,? 

10  D(K)=2.*A 

DO  15  K=1.N.2 
15  0(K)=4.«A 

30  0(K3)=A 

»  CALCULATE  NONHOMOGENFoUS  TERMS 
PHS=1,0 
DO  22  I=1.K2 
print  9 

35  9  FORMAT(lHl) 

DO  999  n  =  l.K4 
DO  35  N=1.K3 

35  NON(N)=RHS»SQRT (PT (N) ) 
call  CONST (I) 

AO  *  CALCULATE  KERNEL  MATRtCFS 

DO  20  N=1»K3 
DO  20  M=1»K3 

F(M,N.I)=FU(I.PT(M) ,PT(N) ) 

20  CONTINUE 

45  call  CHANGE(F.G.D. T ) 

call  L1NEQ(G.8.C»  k3) 

DO  40  L=1.K3 
print  6.PT  (L)  .NONd.) 

6  F0RMAT(5X»F8.4.F15.6) 

50  40  CONT  I  JUE 

LP(1I^1)=N0N(K3) 

DTA(II^1)=P 
999  CONTINUE 

call  LAPINV(DTA.LP) 

55  22  CONTINUE 

END 

-22- 


1  function  SIMP(I,AfP) 

C0MM0N/AUX/H»P»PK1 .PK2*8MU»X,Y 

MXYZ=2**15 

DEL=0.25*{B-A) 

5  IF(DEL)40»45*50 

45  SlMPs0*0 

■  .  RETURN 

50  CONTINUE 

SA=Z(I*A)*Z(If8) 

10  SBsZ(I*A>2,«DEL) 

SC=Z ( I fA^DFL) ♦! ( I ♦ A*3.*OEL) 

SU  (DEL/3.  )*  (SA>2.<»SR>4.*SC) 

IF(S1 .EQ.0.0)  GO  TO  45 
K=8 

15  35  S0=S8*SC 

DEL=0.5*DEL 

SC=Z(I*A*DEL) 

J=K-1 

DO  5  N=3*J*2 
20  AN=N 

5  SC=SC*Z(l»A»AN»OEL) 

S2s (DEL/3. )*(SA>2.*SR*4,*SC) 

0IF=ABS( (S2-S1)/S1> 

ER=0.01 

?5  IF(DIF-ER)30»25.25 

30  SIMP=S2 

RETURN 
25  K=2*K 

Sl=S2 

30  IF (k-MXYZ) 35.35.40 

40  PRINT  42. I. A. 8 

42  F0RMAT(5X.*  INT.  DOES  NOT  CONVERGE  *.I3.2F9.4) 
PRINT  60.X,Y 
60  rORMAT(2F10.5) 

35  DO  70  Jsl.lO 

DiPsJ 

DIP=DIP/10. 

W=Z(I.DIP) 

PRINT  60. W 

40  70  CONTINUE 

call  EXIT 

end 


symbolic 

reference  MAP 

(R=l) 

ENTRY 

POINTS 

4 

SIMP 

variables  SN 

TYPE 

relocation 

0 

A 

REAL 

F.P. 

260 

AN 

real 

0 

6 

REAL 

F.P. 

4 

BMU 

REAL 

250 

del 

REAL 

262 

OIF 

REAL 

264 

DIP 

REAL 

263 

ER 

REAL 
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1 


5 


10 


SUBROUTINE  CHANGE(r*GfDtI) 

COMMON  K1*K2*K3»K4 

peal  F (A,4,2) ,6(A,4) ♦D(4) 

DO  10  N=1»K3 

DO  10  M=1*K3 

G(M,N)  =F(M*N*I)*0(\') 

10  CONTINUE 

DO  20  N=1*K3 
20  G(N,N)=G(N,N)*1.0 

RETURN 
end 


SYMBOLIC  REFERENCE  MAP  (R=l) 


entry  POINTS 
3  CHANGE 


variables 

SN  TYPE 

RELOCATION 

0 

D 

REAL 

ARRAY 

F.P. 

0 

F 

REAL 

0 

G 

REAL 

ARRAY 

F.P. 

0 

I 

INTEGER 

0 

K1 

INTEGER 

/  / 

1 

K2 

INTEGER 

2 

K3 

INTEGER 

/  / 

3 

K4 

INTEGER 

53 

M 

INTEGER 

52 

N 

integer 

statement  labels 

0 

10 

0 

20 

LOOPS 

label 

INDEX 

EROM-TO 

length 

PROPERTIES 

17 

10 

o  N 

4  7 

170 

NOT 

INNER 

30 

10 

M 

5  7 

OB 

INSTACK 

43 

20 

N 

0  9 

aB 

INSTACK 

common 

BLOCKS 

LENGTH 

/  / 

4 

STATISTICS 

program  length  65B  53 

SCM  blank  common  length  4B  4 

47000B  SCM  USED 
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1 


5 


10 


15 


?0 


?5 


30 


35 


40 


SUBROUTINE  LINEO ( A . B » T ^N) 
PEAL  A(N*N) »B(N) ♦T(N) 

DO  5  I=2*N 

5  A(I,1)=A(I,1)/A(1,1 ) 

DO  10  K=2*N 
M=K-1 

DO  15  I=1*N 
15  T(I)=A(I*K) 

DO  20  J=1«M 
A(J*K)s:T(J) 

J1=J*1 

DO  20  I=J1»N 
T(I)=T(I)-A(I,J)»A(J,K) 

20  CONTINUE 

A (K,K)sT (K) 

IF(K*EQ.N)  go  TO  10 
M*K  +  1 

DO  25  I=:M*N 
25  A(I.K)=T(I)/A<K,K) 

10  CONTINUE 
*  BACK  SUBSTITUTE 
DO  31  I=1*N 
T(I)=B(I) 

Msl^l 

IF(M.GT.N)  go  TO  31 
DO  30  J=M»N 
B(J)=B(J)-A(J,I)*T(I) 

30  CONTINUE 
31  CONTINUE 
DO  35  I=1»N 
K=N*1-I 

R(K)=T(K)/A(K«K) 

Kl=K-l 

IF(KI.EG.O)  GO  TO  35 
DO  36  Jlsl.Kl 
J=K-J1 

T(J)=T(J)-A(J,K)*B(K) 

36  CONTINUE 
35  CONTINUE 
RETURN 
END 


SYMBOLIC  REFERENCE  MAP  (R=l) 


entry  points 

3  LINED 


variables 

0  A 

SN  TYPE 

REAL 

RELOCATION 
ARRAY  F.P, 

0 

B 

REAL 

172 

I 

INTEGER 

175 

J 

INTEGER 

176 

J1 

INTEGER 

173 

K 

INTEGER 

177 

K1 

INTEGER 

174 

M 

INTEGER 

0 

N 

INTEGER 

F.P. 

0 

T 

REAL 
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1 


FUNCTION  FU(I*A,B) 
C0MM0N/AUX/H»P»PK1 tPKZtRMUf X,Y 
X=A 
YsB 

5  IF(A»B)5*10*5 

10  FUsO.O 

RETURN 

5  SUMsSIMPdtO. 0*5.0) 

ER=0.01 

10  DEL  *5.0 

20  UP*OEL+5.0 

ADOL=SIMP(I,OEL*UP) 

DEL  *UP 

TEST=ABS(AODL/SUM) 

15  SUM=SUM*AODL 

IF(TEST-ER) 15*20*20 
15  FUsSQRT()<»Y)*SUM 
RETURN 
END 


SYMroLIC  reference  MAP  (Psl) 

ENTRY  POINTS 
4  FU 


VARIABLES  SN  TYPE  RELOCATION 


0 

A 

REAL 

E.P. 

62 

addl 

real 

0 

B 

REAL 

i^.P. 

4 

BMU 

real 

60 

del 

REAL 

57 

ER 

REAL 

55 

EU 

REAL 

0 

H 

REAL 

0 

I 

INTEGER 

E.P. 

1 

P 

real 

2 

PKl 

REAL 

AUX 

3 

PK2 

REAL 

56 

SUM 

REAL 

63 

TEST 

REAL 

61 

UP 

REAL 

5 

X 

real 

6 

Y 

REAL 

AUX 

externals 

TYPE 

ARGS 

SIMP 

REAL 

3 

SORT 

REAL 

INLINE  FUNCTIONS  TYPE  ARCS 

ABS  REAL  1  INTRIN 

STATEMENT  LABELS 

14  5  0  10  INACTIVE 

22  20 

common  BLOCKS  LENGTH 
AUX  7 

STATISTICS 

PROGRAM  length  64B  52 

SCM  LABELED  COMMON  LENGTH  7B  7 

47000B  SCM  USED 
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1  function  8FSJ0(A) 

IF(A-3.)5*5*10 
5  B=A*A/9. 

W*l. -2. 2499997*8 
5  Z=B*B 

W=W*1.2656208*Z 

z*z*e 

WsW-.3163866*Z 

Z=Z*B 

10  W=W+.0444479*Z 

Z=Z»8 

W=W-.0039444*Z 

7=2*8 

BESJO=W*.00021*Z 
15  RETURN 

10  8=3. /A 

W=. 79788456-. 00 OOOn77*B 
V=A-. 785398 16-. 041 66397*8 
Z=B*B 

?0  W=W-.0055274*Z 

V=V-.00003954»Z 

Z=Z*B 

W=W-. 00009512*7 
V=V^.00262573«Z 
25  Z=7*R 

W=W>.00137237*Z 
V=V-. 00054125*2 
Z=Z*B 

W=W-. 00072805*2 
30  V=V-. 00029333*2 

Z=Z*B 

W=W^. 00014476*7 
V=V4. 00013558*2 
BES JO=i</SQRT  ( A )  *CO':  <  V ) 

35  RETURN 

END 


SYMBOLIC  REFERENCE  MAP  (R=l) 

entry  points 

4  BESJO 


variables  SN  type  RELOCATION 


.  0  A 

REAL 

F.P. 

114 

8 

real 

113  BESJO 

PEAL 

117 

V 

REAL 

115  W 

PEAL 

116 

z 

REAL 

EXTERNALS 

TYPE 

ARGS 

COS 

'  REAL' 

1  LIBRARY 

SORT 

real 

statement  LABELS 

0  5  INACTIVE  10 
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SUBROUTINE  CONST ( I ) 

COMMON/AUX/HfP»PKl ,PKH»BMU»X,Y 

PK1=0.3 
PK2=0,3 
BMUxBO.O 
Hsl.O 
read  2»P 
’  FORMAT (FIO. 5) 

hh=i ./h 

PRINT  1»BMU*PK1»PK?»HH,p 

F0RMAT(/////5X*»  mU?/MU1  =*F6.2**  NUl  =»F4,2*»  NU2 
lA/H  =:*F4,2.*  C21/oa  =*F4,2//) 

RETURN 

END 


symbolic  reference  map  (R=1) 


ENTRY  POINTS 
3  CONST 


variables 

SN  TYPE 

relocation 

4  BMU 

REAL 

AIJX 

55  HH 

REAL 

1  P 

REAL 

AUX 

3  PK2 

REAL 

AUX 

G  Y 

REAL 

AUX 

FILE  names 

MODE 

input 

FMT 

output 

0  H 
0  I 
2  PKl 
5  X 


REAL 

INTEGER 

REAL 

REAL 


STATEMENT  LABELS 
37  1  FMT 


25  2 


common  blocks  length 

AUX  7 


statistics 

program  length 
SCM  LABELED  COMMON  LENGTH 
47000B  SCM  USED 
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I 


1  function  Z(I*S) 

COMMON/AUX/H»P»PK1,PK2,9MU»X,Y 
COMPLEX  DA»DL1 *OL2.SA»SB*SC*SD 
COMPLEX  GA*G8»CA*CB»CC,F,G 
5  PI=3. 1415926 

pp=p*p 

P6s2./PP/(1.^BMU) 

AAs2,*(1.-pk1)/(1.-2.*pk1) 

ABs2.* ( 1 •-PK2) / ( 1 .-2.*PK2) 

10  PAs2,/PP/(AA*BMU*AO) 

P0s2.*H*H/PI/Pl/ (AA*BMU*AB)/PP 
PA=(l.*0MU)/{AA4'eMii*AB) 

BB=1,-BA 

BC=(l.*7.*BMU)/4,/(l .♦BMU) 

15  BD=Pl*PI/2./H/H 

ALPsl ./4,/BA/BB 
DLP=BC/4./PB 

DD=(  (ALP^OLP)*PO*l  ,)*»?-4,*ALP*»DLP*P0»(P0*1  .) 
6*CMPLX (DD»0.0) 

20  OAsCSQPT(G) 

DLl=BO/OLP»( (ALP*nLP)*PO*l.*DA) 

0L2sBD/0LP« ( (ALP*0LP) *P0^1 .-OA) 

SC=S*S*DL1 

SD=5*S*0L2 

25  GAsCSQRKSC) 

GB=CS0RT(SD) 

GC=SQRT(S*S*PG) 

SAs:(PA-DL2)/(0L1-0L2) 

SB= (PA-OLl ) / (PA-DL5) 

10  CA=SA/PG/BB 

CB=2.* (S*S*PG/2. ) »«2/GA*  < I .-GA/GB*SB) 

CC=2.*S*S*GC/SA 

F=CA*(C0-CC) 

Q=REAL (F) 

35  0A=AIMA6(F) 

IF(OA-0.0)5»10»5 

10  Z-(Q-S)*BESJO(S*X)i^BFSJO(S»Y) 

RETURN 

5  PRINT  9*P»S*F 

UO  9  FORMAT(4F10.5) 

CALL  EXIT 

end 


SYMBOLIC 

REFERENCE  MAP 

(R  =  l ) 

entry 

POINTS 

4 

Z 

VARIABLES  SN 

TYPE 

RELOCATION 

276 

AA 

REAL 

277 

A8 

real 

306 

ALP 

REAL 

302 

BA 

REAL 

303 

BB 

REAL 

304 

8C 

REAL 

305 

BD 

REAL 

4 

BMU 

REAL 
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1  SUBROUTINE  LAPINV (GLAM , PHI ) 

C  THIS  PROGRAM  EVALUATES  THE  COEFFICIENTS  FOR  SERIES 

C  OF  JACOBI  polynomials  WHICH  REPRESENTS  A  LAPLACE 

C  inversion  INTEGRAL 

5  real  MUL 

dimension  a (50) ♦GLAM(50) *PHI (50) .C(4»50) 
dimension  BK(lOl) ♦TT(lOl) 

C0MM0N/2/TI fTF.DTfMN.BKtTT 
READ  lfNN*MN*MM 
10  1  F0RMAT(3I2) 

READ  2»TI*TF*DT 

2  F0RMAT(3F10.5) 

PRINT  99 

99  format (IHl) 

15  call  SPLICE(GLAM*Pm,MM,C) 

PRINT  101 

101  FORMAT (/////5X**  GLAM  PHI  *) 

PRINT  102*(GLAM(I),PhI(I)*Is1»MM) 

102  FORMAT(5XtF10.5*5X.F10.5) 

20  Mll=MM-l 

PRINT  99 
DO  10  I=1*NN 
READ  3*8ET*DEL 

3  FORMAT(2F10.5) 

25  PRINT  98*BET*DEl 

98  FORMAT (/////5X*»BFtA  »*F5.3»*  DELTA  =»F5.3) 

DO  11  L=1*MN 
ALsL 

S=1./(AL^BET)/0EL 

30  call  SPLINE{GLAM,pv4i.MM,C*S«6) 

F=G*S 

IF(AL-2.)81 ♦82»83 
81  A(l)=(l.*BET)*DEL»f 
GO  TO  11 

35  82  A{2)  =  (  (2.^RET)*0EL<^F-A(1))*(3.*BET) 

GO  to  11 
83  CONTINUE 
TOPrl . 

L1=L-1 

40  AL1=L1 

DO  12  J*ltLl 
AJ=J 

TOP=AJ«TOP 
12  CONTINUE 
45  L2=r2*L-l 

BOT=l. 

DO  13  J=L»L2 
AJ=J 

eOT=(AJ^BET)*BOT 
50  13  CONTINUE 

MUL=B0T/T0P 
SUM=0.0 
DO  14  N=1 tLl 
ANsN 

55  IF(AN-2.)85»86»87 

85  TOO=l. 

GO  TO  88 
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60 


65 


70 


75 


80 


86  TOOsALl 
GO'  TO  88 

87  CONTINUE 
TOO=l, 

lCH=Ll-(N-2) 

DO  15  J=ICH*L1 
AJsJ 

TOD=AJ«TOO 

15  CONTINUE 

88  CONTINUE 

800=1. 

JA=L1*N 

DO  16  J=L»JA 
AJ=J 

BOO=BOD»(AJ*0ET) 

16  CONTINUE 
C0=T00/800 
SUM=SUM^CO*A (N) 

14  CONTINUE 

A(L)aMUL*(DEL*F-SUM) 

11  CONTINUE 

CALL  JACSER(OEL.A.PET) 
10  CONTINUE 
999  CONTINUE 
RETURN 

end 


symbolic  reference  map  (R=1) 


entry  POINTS 
3  LAPINV 


variables 

SN  TYPE 

relocation 

377 

A 

REAL 

ARRAY 

364 

AJ 

real 

3S4 

AL 

REAL 

362 

ali 

real 

371 

AN 

REAL 

351 

BET 

real 

4 

BK 

REAL 

ARRAY  ? 

374 

BOD 

real 

366 

BOT 

REAL 

461 

C 

real 

376 

CO 

REAL 

352 

DEL 

real 

? 

DT 

REAL 

2 

357 

F 

REAL 

3S6 

G 

REAL 

0 

GLAM 

real 

347  ■ 

I 

INTEGER 

373 

ICH 

INTEGER 

363 

J 

INTEGER 

375 

JA 

INTEGER 

353 

L 

INTEGER 

361 

LI 

integer 

365 

L2 

INTEGER 

346 

MM 

INTEGER 

3 

MN 

INTEGER 

2 

344 

MUL 

real 

350 

Mil 

INTEGER 

370 

N 

INTEGER 

34P 

NN 

INTEGER 

0 

PHI 

REAL 

355 

S 

REAL 

367 

SUM 

real 

1 

TF 

REAL 

2 

0 

ti 

REAL 

37? 

151 

TOO 

TT 

REAL 

REAL 

ARRAY  2 

360 

TOP 

real 
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SUBROUTINE  JACSER (n,C.B) 

DIMENSION  C(50) *Sr(50) fPCSO) 

DIMENSION  BK(lOl) .TT(lOl) 

C0MM0N/2/TI ♦TF»DT,MN*0KtTT 

TT<l)sO.O 

BK(1)=0.0 

LM=1 

T=TI 

12  T=T*DT 

Xs2.«EXP(-D*T)-l. 

call  JACOBI (HN*X,B.P) 

SF(l)aC<l)*P(l) 

DO  10  L=2*MN 

LUL-1 

AL=L 

SF(L)sSF(L1 ) ♦C(L)*D(L) 

10  CONTINUE 
lm=lm^i 
BK(LM)=SF(5) 

TT (LM)aT 

IF(T.LE.TF)  go  to  12 
PRINT  97 

97  F0RMAT{/////5Xt*  t  K  T  K 

IT  K  *) 

DO  31  MYsl.25 

MAaMY^l 

MB=MA^25 

MC=mB»25 

mD=mC^25 

PRINT  96*TT (MA) ,BK (MA) ♦TT (MB) tBK (MB) »TT (MC) *RK (MC) *7 
96  F0RMAT(5X^F5.2»3X.c-7.5^3X»F5,2,3Xtr7.5*3X,r5.?,3x,r7 
1F7.5) 

31  CONTINUE 
RETURN 
END 


SYMBOLIC  REFERENCE  MAP  (9=1) 

ENTRY  POINTS 
3  JACSFR 


variables 

SN  TYPE 

RELOCATION 

151 

AL 

PEAL 

0 

B 

real 

4 

BK 

REAL 

ARRAY 

2 

0 

C 

real 

0 

D 

REAL 

F.P. 

2 

DT 

real 

147 

L 

INTEGER 

144 

LM 

integer 

150 

LI 

INTEGER 

153 

MA 

integer 

154 

MB 

INTEGER 

155 

MC 

Integer 

156 

MD 

INTEGER 

3 

MN 

integer 

152 

MY 

INTEGER 

241 

P 

real 

157 

SF 

REAL 

ARRAY 

145 

T 

real 

1 

TF 

REAL 

2 

0 

TI 

real 

151 

TT 

REAL 

ARRAY 

2 

146 

X 

real 
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1  SUBROUTINE  JACOBI  (fg*X  *8*PB) 

C  THIS  PROGRAM  CALCULATES  JACOBI  POLYNOMIALS  OF  ORDER 

C  K-1  WITH  AR6  X  AND  PARAMETER  8  6T  -1 

DIMENSION  PB(N) 

5  AN=N 

IF<aN-2.) 1»2*3 

1  PB(1)*1. 

RETURN 

2  PB(l)rl. 

10  PB(2)sX-B*{1,-X)/2. 

RETURN 

3  BSO=B*B 
BONEsB*!. 

PB(1)=1. 

15  PB(2)*X-B*(1.-X)/2. 

DO  4  K=3*N 
AK=K 

AKUAK-1  . 

AK2SAK-2* 

20  Kl=K-l 

K2=K-2 

C01s( (2.*AKl)^e)»X 

COl=(  {2.*AK2)*B)»C'^1 

COl=( (2.*AK2)*B0NE)*{C01-BSQ) 

?5  C02=2.*AK2»(AK2*8)'><  {2.*AKn*B) 

C0=2.*AK1* (AKl^R) (2.*AK2) >9) 

4  PB(K)=(C01«Pe(Kl)-r02*P8(K2) ) /CO 
RETURN 

END 


SYMBOLIC  reference  MAP  (R=l) 

entry  POINTS 
1  JACOBI 


VARIABLES  SN 

TYPE 

relocation 

10^ 

AK 

REAL 

107 

AKI 

REAL 

1 1  0 

AK2 

REAL 

102 

AN 

real 

0 

B 

REAL 

F.P. 

104 

BONE 

real 

103 

0SO 

PEAL 

115 

CO 

real 

113 

COl 

REAL 

114 

C02 

real 

105 

K 

INTEGER 

111 

K1 

INTEGER 

112 

K2 

INTEGER 

0 

N 

INTEGER 

0 

PB 

REAL 

ARRAY 

F.P. 

0 

X 

real 

STATEMENT  LABELS 

0 

1 

INACTIVE 

24 

2 

0 

4 

LOOPS 

LABEL 

INDEX  FROM-TO 

length 

PROPERTIES 

47 

4 

K 

16  27 

2BB 

OPT 
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SYMBOL!' 
entry  points 

3  SPLINE 


SUBROUTINE  SPL INE ( X ♦ Y *M ,C* X INT ♦ YINT) 
dimension  X(50) .Y(S0) tC(4*50) 

IF(XINT-X(l))l»10,ll 

10  YINT=Y(1) 

RETURN 

11  CONTINUE 
lF(X(M)-XINT)l»12.n 

12  YlNTsY(M) 

RETURN 

13  CONTINUE 
•  KaM/2 

N=M 

2  CONTINUE 
IF(X(K)-XINT)3»14»'= 

14  YINT=Y(K) 

RETURN 

3  CONTINUE 

IF(XINT-X(K>1) )4*ic,7 

15  YlNTsY(K^l) 

RETURN 

4  CONTINUE 

YINT=(X<K*1 )-XlNT)*(C(l .K)*(X (<♦! )-XINT)**2^C(3*K) ) 
YINT=Y1NT^ (XINT-X(k) ) * < C ( 2 ♦ K ) * (X I NT-X (K ) ) **24C (4»K ) ) 
RETURN 

5  CONTINUE 

IF(X (K-1>-XINT)6»1a»17 

6  K=K-1 
60  TO  4 

■  16  YINT=Y (K-l ) 

RETURN 

17  N=K 
K=K/2 
GO  TO  2 

7  LL=K 
K*(N^K)/2 

8  CONTINUE 

IF(X(K)-XINT)3»14,i8 

18  CONTINUE 

IF{X(K-1)-XINT)6«16»19 

19  N=K 
K=(LL^K)/2 
60  TO  8 

1  print  101 

101  FORMAT («  OUT  OF  RAmGE  FOR  INTERPOLATION  *) 

STOP 

END 


i  REFERENCE  map  (R=1) 
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1  SUBROUTINE  SPLICE ( X ♦ Y*M,C) 

DIMENSION  X(50) ♦Y{S0) *0(50) ♦P(50) »E(50) ♦C(4»50) 
DIMENSION  A(50*3) ♦P(50) tZ(50) 

MM=M-1 

5  DO  2  Ksl«MM 

0(K)*X<K^1)-X(K) 

P(K)sO(K)/6. 

2  E<K)=(Y(K*1)-Y(K))/0(K) 

DO  3  Ks2*MM 

10  3  8<K)=E(K)-E(K-1) 

A(l*2)s-1.-D{l)/D(?) 

A(l,3)=0(l)/D(2) 

A(2*3)sP(2)-P(l)*A(l,3) 

A(2*2)=2,«(P(1)>P(P) )-P(l)*A{l*2) 

15  A(2»3)sA(2»3)/A(2,?) 

B<2)=B(2)/A(2*2) 

DO  4  Ks3*MM 

A(K»2)s2.«(P(K-l) ♦P(K) )-P(K-l )»A(K-1*3) 
B(K)sB(K)-P(K-l)*B(K-l) 

20  A(K,3)=P(K)/A(K,2) 

4  8<K)=8(K)/A(K.2) 

Q=0(M-2)/0(M-1) 

A(M*l)»l.*Q*A{M-2,i) 

A(M*2)=-0-A<M»l )*A (M-l *3) 

25  B(M)sB(M-2)-A(M,1)«B(M-1) 

Z(M)=8(M)/A(M»2) 

MN=M-2 

DO  6  I=1»MN 

KsM-I 

30  6  Z(K)=B(K)-A(K»3)«7(K^1) 

Z(l)s-A{l»2)»Z(?)-A  (I 
DO  7  K=1*MM 
0=1 ./(5.*D{K) ) 

C(1  *K)=Z(K)*Q 

35  C(2,K)=Z<K*1)»0 

C(3fK)=Y(K)/D(K)-Z(K)*P(K) 

7  C(4,K)=Y(K*l)/0{K)-Z(K4l)*P(K) 

RETURN 

END 


SYMBOLIC  reference  MAP  (R=l) 

entry  POINTS 
3  SPLICE 


VARIABLES  SN  TYPE  RELOCATION 


373 

A 

REAL 

ARRAY 

621 

8 

REAL 

0 

C 

REAL 

ARRAY 

F.P. 

145 

0 

real 

3]  1 

E 

REAL 

ARRAY 

144 

I 

INTEGER 

141 

K 

INTEGER 

0 

M 

INTEGER 

140 

MM 

INTEGER 

143 

MN 

INTEGER 

227 

P 

REAL 

ARRAY 

142 

0 

real 

0 

X 

PEAL 

ARRAY 

F.P. 

0 

Y 

real 
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